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Abstract 

The ENUF method, i.e., Ewald summation based on the Non-Uniform FFT tech- 
nique (NFFT), is implemented in Dissipative Particle Dynamics (DPD) simulation 
scheme to fast and accurately calculate the electrostatic interactions at mesoscopic 
level. In a simple model electrolyte system, the suitable ENUF-DPD parameters, 
including the convergence parameter a, the NFFT approximation parameter 
and the cut-offs for real and reciprocal space contributions, are carefully deter- 
mined. With these optimized parameters, the ENUF-DPD method shows excellent 
efficiency and scales as 0{N \ogN). The ENUF-DPD method is further validated 
by investigating the effects of charge fraction of polyelectrolyte, ionic strength and 
counterion valency of added salts on polyelectrolyte conformations. The simulations 
in this paper, together with a separately published work of dendrimer-membrane 
complexes, show that the ENUF-DPD method is very robust and can be used to 
study charged complex systems at mesoscopic level. 
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1 Introduction 



The Dissipative Particle Dynamics (DPD) simulation method pip] was origi- 
nally developed to cover much longer length and time scales than in conven- 
tional atomistic Molecular Dynamics (MD) simulations using a mesoscopic 
description of the simulated system [3] due to soft forces acting between 
large particles made of clusters of atoms. Since its introduction about two 
decades ago, several improvements and generalizations |l|^|5|7f5|^lU|ll|12|13] 
have been proposed, making DPD method one of the mostly used coarse- 
grained approach in soft matter simulations. Examples of such studies are 
microphase separation of block copolymers [I^|T5] . polymeric surfactants in 
solution [T6|T7] . colloidal suspensions [18], and the structural and rheological 
behavior of biological membranes [T^PP] . For many of these complex systems 
mentioned above, the electrostatic interactions play a vital role behind the 
key phenomena and dynamical processes. The inclusion of electrostatic inter- 
actions in DPD simulations is important to capture the long-range interactions 
at mesoscopic level of material description for many systems such as the con- 
formational properties of polyelectrolyte brushes pT|22|23] and the formation 
of membrane-DNA complexes in biological systems j24|l25| . 

However, when electrostatic interactions are incorporated in DPD method, 
the mesoscopic DPD particles carrying opposite charges tend to form artifi- 
cial aggregates due to the soft nature of the conservative interactions (forces) 
in DPD simulations. In order to avoid a collapse of oppositely charged parti- 
cles onto each other, Groot [26] smeared out the local point charge into grids 
around each DPD particle. He adopted a variant of particle-particle particle- 
mesh (PPPM) approach, which was originally introduced for systems with 
electrostatic heterogeneities [27], to treat separately the near and far field in- 
teractions between charge distributions. As the charge density distributions in 
the simulated system are affected by hydrodynamic flow [2H], this method pro- 
vides a natural coupling between electrostatics and fluid motion. Simulation 
results [2T|23II25||26] demonstrated that this method is reasonably efficient in 
capturing important features of electrostatic interactions at mesoscopic level. 

In an alternative approach, Gonzalez- Melchor et al. [22] adopted traditional 
Ewald summation method and Slater-type charge density distribution in DPD 
simulations. This method allows a standard approach to calculate electrostatic 
energy and force of charge density distributions. Although the inclusion of 
charge density distribution does not directly increase the computational cost 
in the Ewald summation itself, this method becomes computationally more 
demanding than the one adopted by Groot [26] . 

In order to improve the efficiency and accuracy in treating electrostatic inter- 
actions, we suggest here an alternative approach that allows fast calculation 
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of electrostatic energy and force in DPD simulations. The ENUF method, an 
abbreviation for the Ewald summation using Non-Uniform FFT (Fast Fourier 
Transform) (NFFT) technique, was recently suggested in our group PUIIST] and 
showed excellent computational efficiency in atomistic MD simulations. In our 
current work, we implement the ENUF in the DPD method. The ENUF-DPD 
method is initially applied on simple electrolyte system as a typical model 
case to optimize the ENUF-DPD parameters and investigate corresponding 
computational complexity. With suitable parameters, we adopt the ENUF- 
DPD method to study the dependence of polyelectrolyte conformations on 
ionic strength and counterion valency of added salts for illustration. 

This paper is organized as follows: Sec. [2] contains a brief introduction to the 
DPD method. Detailed algorithm of the ENUF method is given in Sec. |3l 
Sees, m and [5] describe the implementation of ENUF in DPD method and the 
exploration of suitable parameters for ENUF-DPD method in a simple model 
electrolyte system. In Sec. El the ENUF-DPD method is further validated in 
studying polyelectrolyte conformations upon addition of salts with multivalent 
counterions. Finally, main concluding remarks are given in Sec. [3 



2 The DPD Method 



The DPD method, originally introduced by Hoogerbrugge and Koelman in 
1992 [1], is a mesoscopic and particle-based simulation method based on a 
set of pairwise forces. One important conceptual difference between DPD and 
conventional MD is the use of coarse-graining (CG) procedure allowing a map- 
ping of several atoms or molecules from the real atomistic system onto larger 
DPD particles. After scaling up the size of the system, the DPD method can 
capture the hydrodynamic behavior in very large complex systems up to mi- 
crosecond range and beyond. Like in MD simulations, the time evolutions of 
DPD particles are governed by Newton's equations of motion 



— — = Vj , mj— — = ij , (1) 
dt ' dt ' ^ ^ 

where r^, Vj and fj denote the coordinate, velocity, and the total force acting 
on particle i, respectively. The total force, between any pair of DPD particles 
i and j, is normally composed of three different pairwise additive forces: the 
conservative force F^, the dissipative force F^, and the random force F^, 



(2) 



with 
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(3) 
(4) 
(5) 



where r^j = rj — Tj, Vij = \ rij\, fjj = rij/rij, and Vjj = Vj — v^. The parame- 
ters aij, 7, and a determine the strength of the conservative, dissipative, and 
random forces, respectively. 6ij is a randomly fluctuating variable, with zero 
mean and unit variance. 

The pairwise conservative force is written in terms of a weight function a;*"(rjj), 
where u'-^{rij) = 1 — Vij/Rc is chosen for Vij < Rc and uo'^irij) = for Vij > Rc 
such that the conservative force is soft and repulsive. The unit of length Rc is 
related to the volume of DPD particles. In our simulations, we adopt the CG 
scheme [32] with = 4 and p = 4, in which the former parameter means 
4 water molecules being coarse-grained into one DPD particle and the latter 
means there are 4 DPD particles in the volume of R^. With this particular 
scheme, the length unit Rc is given as Rc = 3.107.^^^ = 7.829A. The 
conservative interaction parameters between different types of DPD particles 
are determined by ~ an + 2.05% with an = 78.67/cbT, in which x is the 
Flory-Huggins parameter between different types of DPD particles. 

Unlike the conservative force, two weight functions uj^irij) and uj^ijij) for 
dissipative and random forces, respectively, are coupled together to form a 
thermostat. According to Espanol and Warren [3], the relationship between 
two functions is described as 



This precise relationship between dissipative and random forces is determined 
by the fluctuation-dissipation theorem. We adopt a simple choice of uj^{r) due 
to Groot and Warren [1] 



For polymer and surfactant molecules, the intramolecular interactions between 
bonded particles are described by harmonic springs 




(6) 




Ff = - 




eq 
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where is the spring constant and r^q is the equihbrium bond length. 

In traditional DPD method, charge density distributions are usually adopted 
instead of point charges to avoid the divergence of electrostatic interactions at 
r = 0. In our implementation, a Slater- type charge density distribution with 
the form of 



P.(r) = ^.* (9) 

is adopted, in which Ag is the decay length of charge q. The integration of 
Eq. in] over the whole space gives the total charge q. 

A modified version of velocity- Verlet algorithm |1] is used to integrate the 
equations of motion. For easy numerical handling, we choose the cut-off radius, 
the particle mass, and as the units of the simulating system, i.e., Rc = 
m = ksT = 1. As a consequence, the unit of time r is expressed as r = 



Rcym/kBT = 1. All the related parameters used in our DPD simulations are 
listed in Table [TJ 



3 The ENUF method 



3.1 The Ewald summation method 



Consider a system composed of charged particles, each one carrying the 
partial charge qi at position in a cubic cell with the volume V = . Over- 
all charge neutrality is assumed in the simulations. For simplicity, only the 
simple charge-charge interaction is considered, interactions between dipoles 
and multipoles are omitted in our current scheme. Charges interact with each 
other according to the Coulomb's law, and the total electrostatic energy can 
be written as 



where n = (77,3,, n^, n^), and n^., n^, and are integer numbers. The sum over 
n takes into account the periodic images, and the f symbol indicates that 
the self-interaction terms are omitted when n = 0. The variables eg and 
are the permittivity of vacuum and the dielectric constant of water at room 
temperature, respectively. 
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In the Ewald summation method, the electrostatic energy, as shown in Eq. [101 
is decomposed into real space and reciprocal space contributions [33]. With 
such decomposition, both real and reciprocal space contributions are short- 
range and Eq. [TU] can be rewritten as 



ATxe^er 



EE 



mi 



erfc(^ 



a r,; 



V ^ 

^ n^O 



g-|n|2/4a2 



^(n)5(-n)-^5:gn. 



nL 



N 



n 



(11) 



with 



In 

5(n) = ^gie-*" ""^ and n = — (n^, n^, n^) . (12) 

i=i 

The first, second, and last terms in the bracket on the right-hand side of Eq. [TT] 
correspond to the electrostatic energies from real space, reciprocal space and 
self-interaction parts, respectively, a is the Ewald convergence parameter and 
determines the relative convergence rate between real and reciprocal space 
summations, n is the magnitude of the reciprocal vector n. Choosing suitable 
parameters, the complexity of the Ewald summation method is reduced from 
(9(A^^) to OiN'^l'^^ with considerably accuracy and efficiency [Mj . 

As the number of charged particles in the simulated system grows, it is con- 
venient to combine the calculation of the short-range conservative force with 
the calculation of the real space summations of electrostatic interactions. The 
cut-off for conservative force calculations should be the same as the cut-off for 
real space electrostatic interactions. With such combination and suitable value 
of a, the summation of real space electrostatic energy between two charged 
particles extends no longer than the cut-off distance, and can be expressed as 



^E^R ^ ^_ y y Mierfc(ar) . (13) 

The real space electrostatic force on particle i is the negative of the derivative 
of the potential energy U^'^ respect to its position rj. A common form of the 
real space electrostatic force is described as 



-\E,R _ 




^^erfc(ar) 



2a 



TT 



(14) 
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Thus the real space electrostatic energy and force can be calculated together 
with conservative force, but then the computation of reciprocal space sum- 
mations of electrostatic interactions becomes the more time-consuming part. 
The introduction of the FFT technique reduces the computational complexity 
of Ewald summation to 0{N log N) by treating reciprocal space summations 
with FFT technique [30ll3ll35ll36ll371l38] . 

3.2 Basic features of the FFT technique 

Initially for a finite number of given Fourier coefficients fk^C with k G Im, 
we wish to evaluate the trigonometric polynomial 

fix) = J2 he-'"'"'-^ (15) 
at each of the given nonequispaced points 

X = {x,eD^:j = 0,l,...,N-l} 
which are randomly localized in ci- dimensional domain 

D'' = {x = {xt)t=o,i,...,d-i ■.-l<xt<^}. 

The space of the d-variable function / G D"^ is restricted to the space of 
(i- variable trigonometric polynomials ^e"^'^*'" : k G /m) with degree Mt{t = 
0, 1, . . . , d — 1) in the t-th dimension. The possible frequencies k are collected 
in the multi index set Im with 

h, = {k = {kt)t=o,i,...4~i ^Z^:-^<h<^]. (16) 

The dimension of the function space or the total number of points in the index 
set is Mn = nfjoMt. 

With these in prior definitions, the trigonometric polynomials for the given 
points can be described by 

= /(x,) = (j = 0, 1, . . . , iV - 1) . (17) 

Using the matrix- vector notation, all trigonometric polynomials can be rewrit- 
ten as 
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/ 

where 



= Af, 



(18) 



/ = (/,).=0,l,...,7V-l, (19) 
— l*^ )j=0,l,...,N-l;kelM ^ 

f = {fk)keiM ■ 



In the following implementations, the related matrix-vector products are the 
conjugated form 



f = Af, fj = (20) 



and the transposed form 



N-l 

f = A^f, h=Y. f.e-^^'"-^^ , (21) 

3=0 

in which the matrix A and A^ are the conjugated and transposed complex 
of the matrix A, respectively. With given Fourier coefficients /, the Fourier 
samples / can be transformed with suitable FFT techniques in both directions. 



3.3 Detailed description of the ENUF method 



The ENUF method, which combines traditional Ewald summation method 
with the Non-Uniform fast Fourier transform technique, is a novel and fast 
method to calculate electrostatic interactions. The NFFT [39] is a generaliza- 
tion of the FFT technique [ID]. The basic idea of NFFT is to combine the 
standard FFT and linear combinations of window functions which are well 
localized in both space and frequency domain. The controlled approximations 
using a cut-off in the frequency domain and a limited number of terms in the 
space domain result in an aliasing error and a truncation error, respectively. 
The aliasing error is controlled by the over-sampling factor as , and the trunca- 
tion error is controlled by the number of terms p. As described in Refs. [SUEI], 
NFFT only makes approximations on the reciprocal space part of Ewald sum- 
mation. Hence in the following, we show the detailed procedure to calculate 
the reciprocal space summations of electrostatic energy and force with NFFT 
technique. 
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Recasting the reciprocal space electrostatic energy in terms of Fourier compo- 
nents, the second term of Eq. [11] can be rewritten as 



U''''' = E , „ g(n)g(-n) 

47reoe, ^ |n|2 ^ ' ^ ' 



E -2 5(n)5(-n) . (22) 



A-ne^er 2ttL n 

For fixed vector n, the structure factor S{n) is just a complex number. With 
normalized locations, Xj = rj/L, the structure factor can be expressed with 



TV 

S{n) = . (23) 



It should be noted that the structure factor in Eq. [23]and the transposed FFT 
form in Eq. [21] have similar structures. Suppose that qj is substituting fj, 
the structure factor S{n) is then a three-dimensional case of transposed FFT 
form. By viewing the structure factor S{n) as a trigonometric polynomial /„, 
the reciprocal space electrostatic energy can be rewritten as 



1 1 g-(7rn)2/(«L)2 ^ 



47reoe,. 27rL 



n 



The reciprocal space electrostatic energy is approximated by a linear combina- 
tion of window functions sampled at nonequidistant Mn grids. These grids are 
then input to the transposed FFT, with which one can calculate each compo- 
nent of the structure factor S'(n), and hence the reciprocal space summations 
of electrostatic energy. 

Alternatively, the reciprocal space electrostatic energy in Eq. [22] can be ex- 
pressed by the real and imaginary parts of the structure factor S'(n) as 



1 1 e-(^n)V(«L)2 r ^ . 



1 1 _e-(-")'/("^)M,^ ,27r 

E —o i IE^iCos(— riTi 



|2 



+ lE^»sin(^n-ri)p| 



(25) 



Similarly, the reciprocal space electrostatic force on particle i is the nega- 
tive derivative of the potential energy U^'^ respect to its position and is 
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described as 



Ff'^ = -V,U^'^ 



1 ^ g-(^n)V("i)2 /47rg^ 

L — I 

117^0 



Aireoer 2-kL -n? \ L 



I - sin(^n ■ Ti) Qj cos(^n ■ r^) 
+ cos(^n ■ Ti) Qj sm(^n ■ rj)| 

2^ n ^ <^ sm(— n ■ rj/iel 5(n 



+ cos(^n- r,)/m(5(n))| . (26) 

Since the structure factor S{n) is a complex number, the expression in the 
bracket of Eq. [26] can be written as the imaginary part of a product 



sin(^n-r,)i?e(5(n)) + cos(^n ■ ri)/m(5(n)) = /m|e^*"-''»^(n)|(27) 
Following this expression, Eq. |26] can be expressed with 
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n^O 










I n^O 








I n^O 



1 2g, ^ e-{-n)V{"L)2 ^ 



1 2g, 



(28) 



-(7rn)2/{aL)2 



in which = ^ 'S'(n) with n 7^ 0. Again, Eq. [281 can be considered 

as a three-dimensional case of conjugated FFT form. Assuming n G Im and 
gQ = 0, we can reformulate Eq. [28] into Fourier terms 



1 2g, 



Jmj ^ g„e2'^^"-^4 



Jm(g,) . (29) 
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Thus, one can calculate the reciprocal space electrostatic force on particle i 
using conjugated FFT based on the S'(n) obtained from transposed FFT. 

From the equations shown above, it is clear that by using suitable NFFT 
techniques, the reciprocal space summations of both electrostatic energy and 
force can be calculated from ENUF method. The main features of ENUF 
method are listed in Table |2j 



4 The ENUF-DPD method 



In the traditional DPD formulation, the conservative force F/l- between inter- 
acting particles i and j is a short-range repulsive force, modeling the soft na- 
ture of neutral DPD particles. The electrostatic interactions between charged 
particles are long-range and conservative. In the ENUF-DPD simulations, the 
long-range electrostatic forces and short-range conservative forces are com- 
bined together to determine the thermodynamic behavior of the simulated 
systems [1]. 

When the electrostatic interactions are included in DPD method, the main 
problem is that DPD particles with opposite charges show a tendency to 
collapse onto each other, forming artificial ionic clusters due to the soft na- 
ture of short-range repulsive interactions between DPD particles. In order to 
avoid this, the point charges at the center of DPD particles should be re- 
placed by charge density distributions meshed around particles. Groot [26] 
firstly smeared out the local point charges around regular grids, and then 
adopted a variant of PPPM method to solve the near field and far field equa- 
tions on grids instead of using FFT technique. Gonzalez-Melchor et al. [29] 
directly adopted the Slater-type charge density distribution and traditional 
Ewald summation method to calculate the electrostatic interactions in DPD 
simulations. In our ENUF-DPD method, similar Slater-type charge density 
distribution and Gaussian type window functions in NFFT are adopted, re- 
spectively, to calculate real and reciprocal space contributions of electrostatic 
interactions. 

The detailed procedure of deducing electrostatic energy and force between two 
Slater-type charge density distributions are described as follows. The electro- 
static potential field 0(r) generated by the Slater-type charge density distri- 
bution Pe{f) can be obtained by solving the Poisson's equation 

V20(r) = -— pe(r) , (30) 
In spherical coordinates, the Poisson's equation becomes 
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For the Slater-type charge density distribution in Eq. [9], we define parameter 
c = — for the convenience of integration. Multiply at both sides of Eq. 
and then integrate this equation, we can get 



dr 



eoCr nXl Jo 



2r 



2+~'^ 

C C & ' 



(32) 



By dividing at both sides of Eq. [32], the potential field can be integrated 
analytically to give 
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(33) 



With the definition of c = — j-, we can reformulate Eq. 1551 as 



1 g 



4r 



1 g / (^]^ _|_ ^ — 
Aneoer r ^ Ae ^ 



4r^ 



(34) 



The electrostatic energy between interacting particles i and j is the product 
of the charge of particle i and the potential field generated by particle j at 
position Ti 



(l + ?^)e^). (35) 



Aneoer ^ X^ 



By defining dimensionless parameters r* = r/Rc as the reduced center-to- 
center distance between two charged DPD particles and (3 = Rc/Xg, respec- 
tively, the reduced electrostatic energy between two Slater-type charge density 
distributions is given by 
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The electrostatic force on charged particle i is 



(36) 



Ff'"^"(r.,) = -V.Ug'---(r.,) 




I I] '^e 



11. 



1 _ _ - i> f37) 



With two dimensionless parameters (3 and r*, the magnitude of the reduced 
electrostatic force is 



Comparing Eq. [10] with Eq. [361 we can find that the electrostatic energy be- 
tween two charges is scaled with correction factor i?i = 1 — (1 + /3r*)e~^'^'' 
when the Slater-type charge density distributions are introduced in DPD sim- 
ulations. Similarly, the electrostatic force between two charged particles is 

scaled with correction factor ^2 = 1 — ^1 + 2/3r*(l + /3r*)^e"^^''* in DPD 

simulations. 

In the limit of r*- — )■ 0, the reduced electrostatic energy and force between two 



charge density distributions are described by lim U^'^^'^(r*j) = -^^^'^P> 
and lim F^'^^'^(r* ) = 0, respectively. It is clear that the adoption of Slater- 

type charge density distributions in DPD simulations removes the divergence 
of electrostatic interactions at r*, = 0, which means that both electrostatic en- 
ergy and force between two charged particles are finite quantities. By match- 
ing the electrostatic interactions between two charge density distributions 
at r*j = with previous work [26] gives /3 = 1.125. From the relation of 
/3 = Rc/Xe, we can get Ag = 6.954 A, which is consistent with the electrostatic 
smearing radii used in Gonzalez-Melchor's work [29] . 

Fig. [T] shows the representation of the reduced electrostatic energy and corre- 
sponding force with respect to the distance between two charged DPD parti- 
cles. Both the electrostatic energy and force are calculated using ENUF-DPD 
method and Ewald summation method with Slater-type charge density distri- 
butions and reference parameters. For comparison, we also include the stan- 
dard Coulombic potential and corresponding force, both of which do diverge at 
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r = 0. The electrostatic energy and force calculated using ENUF-DPD method 
are almost the same as those calculated using Ewald summation method with 
reference parameters. The positions of the maximum value of the electrostatic 
energy in both methods are basically the same, but the maximum value calcu- 
lated from ENUF-DPD method is slightly smaller than that from the standard 
Ewald summation method within an acceptable statistical error. Comparing 
with the standard Coulombic potential and corresponding force, we find that 
both ENUF-DPD and Ewald summation methods can give indistinguishable 
energy and force differences at r > 3.0-Rc- Hence, the ENUF method can cap- 
ture the essential character of electrostatic interactions, as well as the Ewald 
summation method, at mesoscopic level [29] . 

Combining the electrostatic force F^'^^^ and the soft repulsive force gives 
the total conservative force F^* between interacting particles i and j in DPD 
simulations. The total conservative force F^*, together with dissipative force 
Ff- and random force F,^, as well as the intramolecular bonding force Ff for 
polymers and surfactants, act on DPD particles and evolve toward equilibrium 
conditions before taking statistical analysis. The details of the update scheme 
for ENUF-DPD method in single integration step are shown in Table El 



5 The choice of ENUF-DPD parameters 

As the number of charged DPD particles in the simulated system grows, 
the calculation of the reciprocal space electrostatic interactions will become 
the most time-consuming part. Using the suitable parameters in ENUF-DPD 
method assures that the time to calculate the real space summations is approx- 
imately the same as the time to calculate the reciprocal space summations, 
thereby reducing the total computational time. Herein we try to explore the 
ENUF-DPD related parameters and get a set of suitable parameters for fol- 
lowing applications. 

The implementation of ENUF-DPD method uses the Ewald convergence pa- 
rameter a, required accuracy 5(<^ 1), and two cut-offs (r^ for real space and 
Uc for reciprocal space). These parameters are correlated with each other with 
the following two conditions 

n^(xL(x N^/^ , 

(39) 

With required 5, it is more convenient to pick a suitable value for ric- Then 
one can determine a and directly from Eq. [391 However, due to the fact 
that Tic should be integer and Tc should be a suitable value for the cell-link 



erfc (arc) ^ 



Tic > 



aL 



TX 



9 9 



14 



list update scheme in DPD simulations, we adopt another procedure to get 
suitable parameters. 

First, we choose suitable 6. It has been demonstrated that 6 = 5.0 x 10^^ is 
enough to keep acceptable accuracy in ENUF method [30]. In DPD simula- 
tions, due to the soft repulsive feature of the conservative force F*^ in Eq. [3l 
we adopt 5 = 1.0 x lO""^ in our ENUF-DPD method. 

Secondly, we determine suitable Vc and a. Gonzalez- Melchor et al. [29] adopted 
1.08-Rc and S.ORc, respectively, as electrostatic smearing radii and real space 
cut-off for the calculation of electrostatic interactions with Ewald summation 
method. In our ENUF-DPD method, as specified in Eqs. [31] and [33 electro- 
static energy jj^^^^^ and force Y^^^^^ are scaled with two correction fac- 
tors, Bi and B2, respectively, both of which are r-dependent. This means that 
the reciprocal space summations of electrostatic energy \je,k,dpd force 
YE,k,dpd g^j,g g^jgQ gQaled with corresponding correction factors. It should be 
noted that what we obtain from FFT is the total influences of other particles 
on particle i. It is difficult to differentiate individual contribution since each 
corresponding correction factor is related to the relative distance between in- 
teracting particles. But if we choose suitable Tc, beyond which two correction 
factors Bi and B2 approximate to 1.0, the total reciprocal space summations of 
electrostatic energy Bi\J^'^'^^^ and force B^F^'^'^^^ can be approximately 
expressed by \jE'^^dpd and y^'^'^^^ , respectively. Such approximation en- 
ables us to adopt directly the FFT results as reciprocal space summations. 
In Fig. [21 we plot two correction factors Bi and B2 with respect to the dis- 
tance r. It is clearly shown that both Bi and B2 approximate to 1.0 when 
r > 3.0-Rc- Hence in our simulations, = 3.0i?c is taken as the cut-off for real 
space summations of electrostatic interactions. With such adoption, both Bi 
and B2 are only applied on real space summations of electrostatic interactions 
within cut-off Tc = 3.0Rc. 

Since the Fourier-based Ewald methods utilize the FFT technique to evalu- 
ate the reciprocal part summations, it is more appropriate to choose suitable 
a, with which we can minimize the total computational time in calculating 
electrostatic interactions. A large value of a means that a small value of Vc is 
used for rapid convergence in real space summations, but the reciprocal space 
calculations will be the more time-consuming part and vice versa. The choice 
of a is system-dependent and related to the trade-offs between accuracy and 
computational speed. Based on Eq. [39] and above determined Vc = S.ORc, we 
deduce that a > 0.12 A . Although the electrostatic energy is invariant to the 
choice of a, the a value indeed affects the total time in calculating electrostatic 
interactions. In order to find a suitable value for a, we carry out our first simu- 
lation to evaluate the Madelung constant (M) of a face-centered cubic (FCC) 
lattice. The crystal structure is composed of 4000 charged particles, half of 
which are cations with net charge +1 and the other half are anions carrying 
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net charge —1. All ions are located on regular lattice of FCC structure. The 
electrostatic energy for the FCC crystal structure is calculated within single 
step and then used to determine Madelung constant M. Simulation details are 
listed in Table El By comparing our calculated M values with the theoretical 
value in Ref. [H], we find that for a wide range of a values the calculated M 
coincides with literature value. The lowest acceptable value, a = 0.20 A , is 
then used in the following simulations to minimize the computational effort. 

Finally, with Eq. [39] and above described parameters, we can choose suitable 
value for ric. For NFFT technique, there are two other parameters, as and 
p, controlling the approximation errors. It has been shown that for a fixed 
over-sampling factor > 1, the error decays exponentially with p [30112] • In 
atomistic MD simulations, it has been shown that cXs = 2 is adequate to take 
enough samples [30|3T] . Hence as = 2 is used in our ENUF-DPD simulations 
to keep a good accuracy. 

Since ENUF approximates only the reciprocal space summations of the stan- 
dard Ewald summation, it is reasonable to expect that both the ENUF and 
the Ewald summation methods should behave in the same way at mesoscopic 
level. We then perform the second set of simulations on bulk electrolyte sys- 
tem to determine parameters, ric and p. The volume for each simulation is 
V = (10-Rc)^ with the total number of charged DPD particles = 4000, in 
which 2000 charged DPD particles represent cations with net charge -|-1 and 
the same number of DPD particles carrying net charge —1 that represent the 
anions. All simulations are equilibrated for 1 x 10^ steps and then another 
2 X 10^ steps to take statistical average for the following analysis. Detailed 
simulation information are listed in Table El 

In Fig. El we present the relative errors of electrostatic energies calculated 
using ENUF-DPD and Ewald summation methods with explored parameters 
ric and p, as well as Ewald summation method with reference parameters. It 
is clear that as Uc increases, the relative errors for electrostatic energy from 
Ewald summation method converge to 1.0 x 10~^, which is the acceptable 
accuracy we set at first. For the ENUF-DPD method, the relative electro- 
static energy errors generally decrease with the increase of p. When p = 1, 
the relative errors from EUNF-DPD method converge to 2 x 10"'^ for large 
ric, indicating that the ENUF-DPD method with p = 1 is too crude due to 
the fact that p is not adequate to provide enough sampling terms in simula- 
tions. By increasing p, the relative errors reduce rapidly with the increase of 
Uc- In Fig. m^a), we show three cases of the electrostatic energy errors calcu- 
lated from ENUF-DPD method with p = 2 and different Uc values, as well as 
those from the Ewald summation method with same simulation parameters. 
It is obvious that for parameters ric > 7, the electrostatic energy errors, calcu- 
lated from both ENUF-DPD and Ewald summation methods, fluctuate within 
1 X 10~^. This means that the electrostatic energy errors from ENUF-DPD 



16 



and Ewald summation method with p = 2 and Uc = 7 are expected to be prac- 
tically negligible in comparison with those calculated from traditional Ewald 
summation method with reference parameters, which gives the most precise 
electrostatic energies. Similar tendencies are also observed in the electrostatic 
force errors calculated from ENUF-DPD and Ewald summation methods with 
the same set of parameters, as shown in Fig. Hl^b). In Fig. |H^c) we show the 
maximum errors in electrostatic forces. It is clear that as p = 2, the increase 
of Uc will enable ENUF-DPD method to be more and more similar to Ewald 
summation method. With larger p values, such as p = 3, we could further 
increase the accuracy of electrostatic interactions in ENUF-DPD method, but 
at the same time the total computational time in treating electrostatic inter- 
actions increases. By compromising the accuracy and computational speed in 
the ENUF-DPD method, we adopt p = 2 and Uc = 7 in following simulations. 

Now we perform the third set of simulations to study the structural properties 
of the electrolyte solution using Ewald summation method with reference pa- 
rameters, ENUF-DPD and Ewald summation methods with above determined 
parameters, respectively. The system consists of = 4000 DPD particles in a 
simulation box of volume V = (10-Rc)^- The solvents mimicking water at room 
temperature are represented by 3736 neutral DPD particles, and 132 particles 
representing ions with net charge -|-1 and the same number of particles with 
net charge —1 representing counterions are also included in the simulations. 
Mapping these quantities to real units, the simulated system corresponds to 
dilute electrolyte solution with a salt concentration about 0.6 M, which is 
consistent with the simulation condition in Refs. Detailed simulation 

information can be found in Table |3l 

The structural properties of neutral solvents and charged particles are de- 
termined by the radial distribution functions (RDFs). The RDFs of neutral 
solvent-solvent, equal sign ions, and unequal sign ions are calculated from 
traditional Ewald summation with reference parameters, the Ewald summa- 
tion and ENUF-DPD method with above determined parameters, as shown 
in Fig. [5l It is clear that the RDFs for the same pair particles calculated from 
three methods show similar tendencies. A general observation is that there is 
no ionic cluster formation at distance close to r = 0. Furthermore, we also find 

that the RDFs between charged particles satisfy g+^{r)g^^/ (r) = gQ^^r), 

where — , and correspond to positive charged, negative charged, and neu- 
tral DPD particles in simulations. This relationship between three RDFs im- 
plies that the structures between charged particles are related to the effective 
electrostatic interparticle potentials, as observed in Ref. 

Another important property of the ENUF-DPD method is its complexity in 
treating electrostatic interactions. Theoretically, the complexity for treating 
electrostatic interactions with FFT related technique is 0{N log N) [31]. For 
NFFT, the computational complexity is 0{MulogMn + log{N/6)), where Mn 
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is the total number of points in the index set, and S is the desired accuracy 
and also a function of p for fixed over-sampling factor as in NFFT [39P2] . 
Combining the definition of Mn and the relationship in Eq. 133 we can get 
Mji (X (X N, hence the theoretical complexity of ENUF-DPD method is 
0{N\ogN + \og{N/6)). 

The scaling behavior of the ENUF-DPD method, as well as the Ewald sum- 
mation method, including its traditional version and the one with suitable 
parameters we have explored above, are estimated by varying the number of 
charged DPD particles in simulations. In each simulation, the total number of 
DPD particles is iV = 32000 with the volume V = {20Rcf. Initially, all DPD 
particles are neutral, and corresponding computational time is taken as the 
benchmark for DPD simulations without electrostatic interactions. Then we 
perform a number of different simulations increasing the number of charged 
particles up to 24000 while keeping the simulation system neutral and the 
total number of DPD particles fixed. The detailed simulation information can 
be found in Table [31 Fig. [6] shows the averaged run time per 10^ steps of 
DPD simulation as function of the number of charged particles. The averaged 
execution time per 10^ steps is the net time in calculating the electrostatic 
energy and force in each DPD simulation. It reveals that the original Ewald 
summation method with reference parameters could generate accurate electro- 
static energy and force, but its computational complexity is 0{N'^). With our 
above determined suitable parameters, the scaling behavior of Ewald summa- 
tion method is reduced to 0{N^^'^) in general. Although the parameters for 
Ewald summation method in our simulations are a little different from the 
parameters used by Gonzalez- Melchor et al. [22] , the computational complex- 
ities are described by similar scaling behavior 0(iV^/^), which has also been 
verified by Ibergay et al. [21]. The ENUF-DPD method scales as C(A^log7V), 
which is in line with the scaling behavior of PPPM method adopted in Groot's 
work [2T1I26] . 

The ENUF-DPD method with above explored parameters shows an excellent 
computational efficiency and good accuracy in treating electrostatic interac- 
tions. In the current study of electrolyte solution and the range of the number 
of ion pairs investigated here, the ENUF-DPD method performs clearly much 
faster than the Ewald summation method, and shows similar 0{N\ogN) scal- 
ing behavior as PPPM method pT|l^ at mesoscopic level. 



6 Conformation of polyelectrolyte in solution 

The above implemented ENUF-DPD method has all capabilities of ordinary 
DPD, but includes applications where electrostatic interactions are essential 
but previously inaccessible. One key example is the polyelectrolyte conforma- 
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tion. Electrostatic interactions between charged particles on polyelectrolyte 
lead to rich conformation of polyelectrolyte qualitatively different from those 
of uncharged polymers [13111] • In this section, we use the ENUF-DPD method 
to study the charge fraction of polyelectrolyte on the conformational behavior 
of single polyelectrolyte molecule. 

Nine charge fractions, defined by / = Nq/Np, in which Nq and Np are, re- 
spectively, the number of charged particles and the total number of particles 
on polyelectrolyte, are considered. In our simulations, Np = 48 is used. Al- 
though this generic polyelectrolyte model is much smaller than the real ones, 
the essential physical effects are still captured [l3|45f46|H7] . Going from a neu- 
tral polymer chain to a fully ionized polyelectrolyte, / takes 0.0, 0.125, 0.25, 
0.375, 0.5, 0.625, 0.75, 0.875 and 1.0, which correspond to Nq = 0, 6, 12, 
18, 24, 30, 36, 42 and 48 charged particles on polyelectrolyte, respectively. 
Each charged particle on polyelectrolyte is characterized by net charge — 1. 
Counterions carrying net charge +1 are added to preserve charge neutrality of 
simulation systems. The equilibrium bond distance between bonded particles 
is set to = 0.7 Rc and the spring constant is taken as = 64.0k bT. The 
conservative interaction parameters between different types of DPD particles 
are obtained through a.^ ^ an + 2.05% with an = TS.GTksT. The x parame- 
ter between polyelectrolyte and solvents is set to 0.87, which is rescaled from 
Refs. [261129] . 

All simulations are performed in a simulation cell with volume V = {30Rc)^. 
The total density is fixed at p = 4, and hence the total number of particles 
in the system is = 108000 in all cases. All simulations are equilibrated in 
5 X 10^ time steps, and then 2.5 x 10^ time step simulations are further per- 
formed to collect statistical data. For a fully ionized polyelectrolyte, additional 
simulations with larger volume V = (40-Rc)^ are also performed. No differences 
beyond statistical uncertainties are found between the two sets of simulations. 
In the following discussion, all simulation results are calculated from systems 
with volume V = {SORcY. Simulation details are listed in Table [31 

The averaged radius of gyration <Rg> of polyelectrolyte, as function of cor- 
responding charge fraction /, are shown in Fig. [71 Experimentally, it is well 
known that a neutral polymer in solution has the smallest radius of gyra- 
tion [^49|50j . When partial groups on polyelectrolyte are ionized, such as 
weakly charged polyelectrolyte in solution with adjustable pH values, the 
<Rg> of polyelectrolyte increases with increasing degree of ionization of poly- 
electrolyte [H]. For fully ionized polyelectrolyte, the value of <Rg> is approxi- 
mately 1.31 times larger than that of neutral polymer, which is consistent with 
molecular simulation results [29] and experimental observations [H]. Typical 
conformations of polyelectrolyte with charge fraction / = 0.0, 0.25, 0.5, 0.75 
and 1.0 are shown in Fig. [HI These simulation results are qualitatively consis- 
tent with experimental observations [ITPS] and theoretical predictions [ISIISU] 
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for weakly charged polyelectrolyte. 



With the increase of charge fraction / on polyelectrolyte, the RDFs between 
charged particles on polyelectrolyte and monovalent counterions are also en- 
hanced, as shown in Fig. [9|^a). In Fig. Mjo), we show the intramolecular pair 
correlation functions between charged particles of polyelectrolyte. A quanti- 
tative measure of the structure is found by analyzing the intramolecular pair 
correlation functions. For neutral polymer and polyelectrolyte with various 
charge fractions, the intramolecular correlations in initial zone r/Rc < 1 are 
dominated by particle-particle repulsions. In the regime r/Rc > 1, two striking 
tendencies are shown in Fig. [9](b). For polyelectrolyte with small charge frac- 
tion, i.e., f < 0.25, we observe a small scaling-like domain and then followed 
by a terminal correlation range. In contrast, polyelectrolyte with large charge 
faction shows a scaling behavior over the entire range. The slope of the fitting 
in Fig. [n](b) is —1.92, which is in good agreement with Groot's results [SB] . 

When salts are added into solution, both ionic strength and valency of mul- 
tivalent counterions of added salts can severely influence the conformational 
properties of polyelectrolyte due to the strong correlations between multi- 
valent counterions and polyelectrolyte. This behavior is usually specified as 
the overcharging phenomenon that occurs in many biological and synthetic 
polyelectrolytes [51]. Herein, the ENUF-DPD method is further adopted to 
investigate the effects of ionic strength and counterion valency of added salts 
on the conformation behavior of fully ionized polyelectrolyte. 

The number of multivalent counterions of added salts, Nc, is determined by 
6 = qNc/Np, where 6 is the ratio between the total charge of multivalent 
counterions of added salts and that of polyelectrolyte, and q is the valency 
of added salts. In our simulations, various 6 values, together with g = 1, 2, 
and 3, are selected to consider the dependence of fully ionized polyelectrolyte 
conformation on the ionic strength and valency of multivalent counterions of 
added salts. All simulations are equilibrated in 5 x 10"^ time steps, and 2.5 x 10^ 
time step simulations are further performed to collect statistical data. Detailed 
simulation information are listed in Table [3l 

Fig. [TU] shows the dependence of <Rg> on 6 and valency of multivalent 
counterions of added salts. In the absence of added salts, polyelectrolyte 
adopts an extended conformation, owing to the electrostatic repulsions be- 
tween charged particles of polyelectrolyte. Upon addition of salts, these re- 
pulsions are screened and hence <Rg> decreases. For monovalent counterions 
(g = 1), <Rg> gradually decreases. By contrast, stronger decreases in <Rg> 
are observed in the cases of added salts with divalent {q = 2) and triva- 
lent counterions (g = 3), occurring at considerably lower ionic strength. This 
reflects the conformational collapse of polyelectrolyte in solution with multiva- 
lent counterions, which has been observed in experiments [ISIIS^ and predicted 
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by coarse-grained molecular dynamic simulations [53]. 

The smallest value of <Rg> occurs near cz, i-e., the (Z:l) salt concentration 
at which the total charge of the Z-valent counterions of added salts neutral- 
izes the bare polyelectrolyte. At the same time, polyelectrolyte shows compact 
conformation. Accordingly, this compact state occurs at a salt concentration 
that decreases with the increase of corresponding counterion valency, which is 
consistent with the two-state model |5l]. In addition, multivalent counterions 
with higher valency are strongly correlated with polyelectrolyte, which can be 
specified by corresponding pair correlation functions. The RDFs among poly- 
electrolyte, monovalent counterions of polyelectrolyte, and multivalent coun- 
terions of added salts at 9 = 1.0, are calculated and shown in Fig. [11] In 
solutions, monovalent counterions of polyelectrolyte and multivalent counte- 
rions of added salts show different condensation abilities on polyelectrolyte. 
For (1:1) salt, two kinds of counterions show similar tendencies due to the 
same amount of net charge on them. With the increase of counterion valency 
of added salts, counterions with different valency show competition in con- 
densating the polyelectrolyte. The peak of RDF between trivalent counterions 
and polyelectrolyte is much higher than that of other counterions, implying 
the strong condensation between trivalent counterions and polyelectrolyte. 
The strong condensation induced by electrostatic correlations decreases the 
osmotic pressure, and hence leads to the collapse of polyelectrolyte [55]. Typi- 
cal conformations of polyelectrolyte, as well as added salt with the counterion 
valency q = 1, 2, and 3, are shown in Fig. [T2J 

A striking effect occurs once the salt concentration is increased beyond cz- 
Polyelectrolyte starts to swell, in close analogy with the redissolution observed 
for multichain aggregates [15] . Comparing with conformation of polyelectrolyte 
in (1:1) salt, which exhibits a slow, monotonic decrease of <Rg> with the 
increase of salt concentrations, the slight swelling behavior of polyelectrolyte 
in the presence of multivalent counterions emphasizes the important role of 
counterion valency [16|52] . 

Concerning the effects of ionic strength and valency of added salts on polyelec- 
trolyte conformation, the decay of correlations between charge particles can 
be specified by the Debye screening length. The addition of salts with multi- 
valent counterions leads to a short Debye screening length [56] , demonstrating 
that the electrostatic interactions between charge particles separated larger 
than specific distance become screened and hence are no longer long-range. It 
is very likely that a finite cut-off for electrostatic interactions, or a screened 
interaction potential between charge particles, can be used in handling elec- 
trostatic interactions in these cases. This topic is beyond the scope of present 
work, deserving a special attention and consideration in future. 

The other detailed application of the ENUF-DPD method is the specific bind- 
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ing structures of dendrimers on amphiphilic membranes [57]. We construct 
mutually consistent coarse-grained models for dendrimers and lipid molecules, 
which can properly describe the conformation of charged dendrimers and the 
surface tension of amphiphilic membranes, respectively. Systematic simula- 
tions are performed and simulation results reveal that the permeability of den- 
drimers across membranes is enhanced upon increasing dendrimer sizes. The 
negative curvature of amphiphilic membrane formed in dendrimer-membrane 
complexes is related to dendrimer concentration. Higher dendrimer concen- 
tration together with the synergistic effect between charged dendrimers can 
also enhance the permeability of dendrimers across amphiphilic membranes. 
Detailed descriptions of this work are shown in Ref. [57] . 



7 Conclusion 



The ENUF method, which combines the traditional Ewald summation method 
with NFFT technique to calculate the electrostatic interactions in MD simu- 
lations, is incorporated in DPD method. The ENUF-DPD method is applied 
on simple model electrolyte systems to explore suitable parameters. With re- 
quired accuracy parameter 6 = 1.0 x 10~^ and cut-off Vc = S.ORc for real space 
summations of electrostatic interactions, we find that the Ewald convergence 
parameter a = 0.20 A can generate accurate Madelung constant for FCC 
lattice structure and keep considerable accuracy. Simulation results reveal that 
the ENUF-DPD method with approximation parameter p = 2 in NFFT and 
cut-off Uc = 7 for reciprocal space summations of electrostatic interactions can 
well describe the electrostatic energy and force, as well as the Ewald summa- 
tion method does. The computational complexity of ENUF-DPD method is 
approximately described as 0{N log N), which shows remarkably better effi- 
ciency than the traditional Ewald summation method with acceptable accu- 
racy in treating long-range electrostatic interactions between charged particles 
at mesoscopic level. 

The ENUF-DPD method is further validated by investigating the influence of 
charge fraction of polyelectrolyte on corresponding conformational properties. 
Meanwhile, the dependence of the conformations of fully ionized polyelec- 
trolyte on ionic strength and valency of added salts are also studied. These 
applications, together with a separately published research work on the forma- 
tion of dendrimer-membrane complexes, show that the ENUF-DPD method is 
very robust and can be used to study charged complex systems at mesoscopic 
level. 
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Table 1 

Summary of the ENUF and DPD parameters used in this work. 

Parameter Value Meaning 

DPD parameters 

Nm 4 Number of water molecules represented by one DPD particle 

p 4 Number of DPD particles in the volume of 

Rc 7.829A Range of soft repulsive interaction in DPD method 

an 78.67 Maximum repulsion parameter between identical DPD particles 

Xij 0.87 Flory-Huggins parameter between polyelectrolyte and solvents 

64.0 Spring constant for bonded interactions of polyelectrolyte 

r^q 0.7 Equilibrium bond length of polyelectrolyte 

7 6.74 Dissipation strength 

a 3.67 Noise amplitude 

A 0.65 Velocity prediction parameter in integration algorithm 

5t 0.02r Integration time step 

m 1 Reduced particle mass 

ksT 1 Energy unit 

ENUF parameters 

5 10~^ Accuracy of electrostatic interactions 

Tc 3.0i?c Cut-off for real space summations of electrostatic calculation 

a 0.20A ^ Ewald convergence parameter 

Us 2 Over-sampling factor 

p 2 Approximation parameter in NFFT 

Uc 7 Cut-off for reciprocal space summations of electrostatic calculation 

Ag 6.954A Charge decay length 

eoEr 1 Dielectric constants 
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Table 2 

The modified version of velo city- Ver let algorithm for the ENUF-DPD method in 
single integration step. 



(0) vO ^ V, + iFf*At + Ff At + F^VAt^ 

(1) vi ^ + (Ff*At + Ff At + Ff VAt) 

(2) Ti ^ ri +v,At 

(3) Calculate Ff {r,}, Ff {r„ v^}, Ff {r,} and Ff{r,} 

(4) Calculate Ff '^^^ and U^'^^^ 

(i) Calculate Ff'^ and U^'-^ 

(ii) Calculate Ff and U^'^ 

5(n) = /„ ^ Eili (using transposed NFFT) 

gn ^ S{n) 

Ei ^ EnG/M gnC^"*" ''" (using conjugated NFFT) 
Ff ^ /^(g.) 

(iii) Calculate U^''^^ (Self-interaction energy) 

(iv) U^.i)^'^ ^ 5i(U^'« + U^'^ - 

(v) Ff'^^^^i?2(Ff'^ + FP) 

(5) Ff*^Ff + Ff + FP^^ 

(6) ^ V, + (fP *At + Ff At + Ff VAt' 
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Table 3 

Simulation details in this work. The 1*** set of simulations is conducted to calculate 
Madelung constant M and determine Ewald convergence parameter a. All charged 
particles are generated on FCC lattice. The 2""^ set of simulations is carried out 
to determine two ENUF parameters, p and ric- The 3'"'^ and 4*'^ simulations are 
used to calculate the RDFs between neutral and charged DPD particles, and to 
investigate the scaling behavior of the ENUF-DPD method, respectively. The 5*'* 
and 6*^^ simulations are performed to study the effects of charge fraction of poly- 
electrolyte, ionic strength and counterion valency of added salts on polyelectrolyte 
conformations, respectively. Representations: N'^ , total number of DPD particles 
in simulation; N~^, number of cations; N~ , number of anions; A'^", number of neu- 
tral particles; /, charge fraction of polyelectrolyte; Np, total number of particles on 
polyelectrolyte; 9, charge ratio between multivalent counterions of added salts and 
that of polyelectrolyte; q, counterion valency of added salts. 
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Fig. 1. Electrostatic potential and force between two charged DPD particles are 
calculated from ENUF and Ewald summation with reference parameters. For com- 
parison, the standard Coulombic potential and force, both of which diverge at r = 0, 
are also included. Both the potential and the force expressions are plotted for two 
equal sign charge density distributions. 

Fig. 2. Two correction factors, Bi and B2, are plotted with respect to the distance 
r. The dotted horizontal and vertical lines imply the value of 1.0 and the cut-off 
Tc = S.ORc for real space summations of electrostatic interactions, respectively. 

Fig. 3. The relative errors of electrostatic energy calculated using the ENUF-DPD 
and Ewald summation with varied p and ric- The relative errors are defined as 

I T T T T^^f f 

I "^ej I , where JJ'^^J is the total electrostatic energy calculated from Ewald sum- 
mation method with reference parameters, U is the total electrostatic energy calcu- 
lated via either the ENUF-DPD or Ewald summation with corresponding explored 
parameters. 



Fig. 4. The errors in electrostatic energy (AU), force (AF) and corresponding maxi- 
mum values (max|AF|) calculated using ENUF-DPD (red squares) and Ewald sum- 
mation (green circles) with p = 2 and various Uc, are compared with those calculated 
from traditional Ewald summation with reference parameters, (a) AU = u— 
where U^e/ is the total electrostatic energy calculated from Ewald summation with 
reference parameters, U is the total electrostatic energy calculated via either EN- 
UF-DPD or Ewald summation with determined parameters, (b) AF = 
where Y^^'^'^f is the averaged electrostatic force on DPD particles calculated from 
Ewald summation with reference parameters, is the average electrostatic force 
on charged particles calculated via either ENUF-DPD or Ewald summation with 

—Y^'^'^f E ref 

determined parameters, (c) max|AF| = max| — ^— jj^^^j — |i=i 2 ■■■ N, where F^ ' ■'is 

the electrostatic force on particle i calculated from Ewald summation with refer- 
ence parameters, is the electrostatic force on particle i calculated via either 
ENUF-DPD or Ewald summation with determined parameters. 

Fig. 5. The pair correlation functions between different types of DPD particles cal- 
culated from Ewald summation with reference parameters, ENUF-DPD and Ewald 
summation with determined parameters. 
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Fig. 6. The scaling behavior of ENUF-DPD method with determined parameters in 
our simulations. The averaged time per 10'^ steps is calculated form DPD simulations 
carried out over 2 x 10^ steps. For comparison, we also show the scaling behavior 
of Ewald summation method, including its traditional version with reference pa- 
rameters and the one with parameters the same as those used in our ENUF-DPD 
method. 



Fig. 7. The averaged radius of gyration of polyelectrolyte as function of charge 
fraction / going from neutral polymer chain to fully ionized polyelectrolyte. 



Fig. 8. Typical conformations of polyelectrolyte with different charge fraction /. 
Solvent particles are not shown for clarity. The red and cyan spheres in all images 
indicate the charged and neutral particles of polyelectrolyte, respectively. In addi- 
tion, counterions are represented by yellow spheres, (a) / = 0.0, (b) / = 0.25, (c) 
/ = 0.5, (d) / = 0.75, and (e) / = 1.0. 



Fig. 9. (a) The intermolecular RDFs between charged particles of polyelectrolyte 
and corresponding monovalent counterions. (b) The intramolecular RDFs between 
charge particles of polyelectrolyte. 



Fig. 10. The dependence of <Rg> of polyelectrolyte on 9 and valency of multivalent 
counterions of added salts. The dash curve gives the <Rg> of neutral polymer in 
solution without salts. 



Fig. 11. The pair correlation functions between polyelectrolyte and counterions, 
including monovalent counterions of polyelectrolyte and multivalent counterions of 
added salts. 



Fig. 12. Typical conformations of polyelectrolyte with added salts at ^ = 1. Solvent 
particles are not shown for clarity. The red and yellow spheres indicate the charged 
particles of polyelectrolyte and corresponding counterions, respectively. Monovalent, 
divalent and trivalent counterions (cations) of added salts are represented by purple, 
green and blue spheres, respectively. All anions of added salts are presented by 
magenta spheres, (a) (1:1) salt, (b) (2:1) salt, and (c) (3:1) salt. 
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Fig. 10. Wang et al. FIGURE 10 
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Fig. 12. Wang et al. FIGURE 12 
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